Load in Precip data

load('data/DailyP.RData')

head(daily_p)
## # A tibble: 6 x 8
##   date       station   lon   lat mean_temp min_temp max_temp daily_p
##   <date>     <chr>   <dbl> <dbl>     <dbl>    <dbl>    <dbl>   <dbl>
## 1 2018-10-01 04V     -106.  38.1      55.7     46.4     71.6    0   
## 2 2018-10-01 0CO     -106.  39.8      38.5     33.8     46.4    0.04
## 3 2018-10-01 20V     -106.  40.1      54.7     39.2     69.8    0   
## 4 2018-10-01 2V5     -102.  40.1      50.9     44.6     68      0   
## 5 2018-10-01 2V6     -103.  40.1      48.7     42.8     68      0   
## 6 2018-10-01 33V     -106.  40.8      54.4     37.4     66.2    0

Get Elevation Data

unique_asos <- daily_p %>%
  distinct(lon, lat, station)  %>%
  st_as_sf(., coords = c('lon','lat'), crs = 4326) %>%
  get_elev_point()
## Downloading point elevations:
## Note: Elevation units are in meters
st_read('data/unique_asso_elev.gpkg')
## Reading layer `unique_asso_elev' from data source 
##   `C:\R-4.1.3\ESS580A9\Metals\data\unique_asso_elev.gpkg' using driver `GPKG'
## Simple feature collection with 67 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -108.7593 ymin: 37.1515 xmax: -102.241 ymax: 40.7503
## Geodetic CRS:  WGS 84

Get Monthly P Averages

monthly_p <- daily_p %>%
  mutate(month = month(date)) %>%
  group_by(month, station) %>%
  summarize(monthly_p = sum(daily_p)) %>%
  left_join(unique_asos) 
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"

Look at monthly P

ggplot(monthly_p, aes(x = elevation, y = monthly_p, color = month)) + 
  scale_color_viridis_c() + 
  geom_point()

Getting Monthly Means of means, mins, maxes in Temp.

monthly_t <- daily_p %>%
  mutate(month = month(date)) %>%
  group_by(month, station) %>%
  dplyr::select(-lon, -lat) %>% 
  summarize(across(where(is.numeric), mean)) %>%
  left_join(unique_asos,.)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"

Temp vs Elevation

ggplot(monthly_t, aes(y = mean_temp, x = elevation, color = month)) + 
  geom_point() + 
  scale_color_viridis_c()
## Warning: Removed 31 rows containing missing values (geom_point).

Assignment

Pick a month (summer months are safer)

July_p <- daily_p %>% 
        st_as_sf(., coords = c('lon','lat'), crs = 4326) %>% 
        mutate(month = month(date)) %>% 
        filter(month == 7) %>% 
        group_by(month, station) %>% 
        summarize(precip = sum(daily_p))
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.

Build IDW precip or elevation for state for that month

# Class hint
# get the boundary data of Colorado
CO <- spData::us_states %>% 
      filter(NAME == 'Colorado')

# make a empty grid
CO_2163 <- st_transform(CO, crs=2163)

CO_2163_box <- st_bbox(CO_2163) %>% 
               st_as_stars(dx = 1000) %>% 
               na.omit(.)

# fit the coordinate reference system
July_p <- July_p %>%
          st_transform(., st_crs(CO_2163_box)) %>% 
          na.omit(.)

# interpolation  
interp_July_p <- idw(precip ~ 1, July_p, CO_2163_box) %>% 
                 dplyr::select(1)
## [inverse distance weighted interpolation]

Plot this data

# using tmap
tm_shape(interp_July_p) +
  tm_raster(palette = 'plasma', style = 'cont',
            title="Total precipitation \n(in) in July 2019") +
tm_legend(legend.outside=TRUE)

# using mapview
mapview(July_p, zcol='precip') + 
mapview(interp_July_p, na.col=NA, 
        col.regions = mapviewGetOption("vector.palette"))

Build IDW with elevation for state for that month including elevation as a predictor

Hint! Use get_elev_raster

# get elevation raster and make grid
ras <- get_elev_raster(unique_asos, z = 8) %>% 
       raster::crop(., unique_asos)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
ras.star <- st_as_stars(ras, dx = 1000) 

names(ras.star) <- "elevation"

# join elevation and precipitation data
unique_asos_July <- st_read('data/unique_asso_elev.gpkg')
## Reading layer `unique_asso_elev' from data source 
##   `C:\R-4.1.3\ESS580A9\Metals\data\unique_asso_elev.gpkg' using driver `GPKG'
## Simple feature collection with 67 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -108.7593 ymin: 37.1515 xmax: -102.241 ymax: 40.7503
## Geodetic CRS:  WGS 84
July_pe <- daily_p %>% 
           mutate(month = month(date)) %>% 
           group_by(month, station) %>% 
           filter(month == 7) %>%
           summarise(precip = sum(daily_p)) %>%
           left_join(unique_asos_July,.) %>% 
           na.omit(.)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"
# Interpolation
interp_July_p_elev <- idw(precip ~ elevation, July_pe, ras.star) %>% dplyr::select(1) 
## [ordinary or weighted least squares prediction]

Make a Map of that

## You will need to create a Stars raster that has elevation data. 

tm_shape(interp_July_p_elev) + 
  tm_raster(n = 10, palette = "-plasma", midpoint = TRUE,
            title="Total precipitation \n(in) in July 2019") +
tm_legend(legend.outside=TRUE)
## stars object downsampled to 1346 by 743 cells. See tm_shape manual (argument raster.downsample)

Compare both maps to PRISM approach for your month

How close do our simple approaches come to reproducing prism maps?

https://www.prism.oregonstate.edu/recent/monthly.php